function F = solve_eqm_inner(param, const, integ)
%% Initial Equilibrium

    % normalization
    P_s     = param.P_s;
    e       = param.e;     
    X       = param.X;

    % read in parameters
    kappa   = param.kappa;
    alpha_m = param.alpha_m; 
    alpha_s = param.alpha_s;
    beta_v  = param.beta_v; 
    beta_m  = param.beta_m;    
    sigma   = param.sigma; 
    sigma_s = param.sigma_s;
    f_vH    = param.f_vH;
    f_vF    = param.f_vF;
    mu_E    = param.mu_E;
    sigma_E = param.sigma_E;   
    w       = param.w;
    Q       = param.Q;
    Nf      = param.Nf;
    
    % read in constants
    gamma   = const.gamma;
    k0      = const.k0;
    C_m     = const.C_m;
    C_v0    = const.C_v0;
    C_x0    = const.C_x0;   
    phi_v   = const.phi_v;
    phi_y   = const.phi_y;
    phi_s   = const.phi_s;
    D_F     = const.D_F;
    lnzQ    = const.lnzQ;
        
    % read in all the integral from previous round
    E_zm0 = integ.E_zm0;
    E_zm1 = integ.E_zm1;
    E_zv0 = integ.E_zv0;
    E_zv1 = integ.E_zv1;
    E_zx0 = integ.E_zx0;
    E_zx1 = integ.E_zx1;
        
    % adjust for level in case integrals too large
    rescale = max(1,10^floor(log10(max(E_zx1))-1));
    E_zx0 = E_zx0/rescale;
    E_zx1 = E_zx1/rescale;
    E_zv0 = E_zv0/(rescale^(1/beta_v));
    E_zv1 = E_zv1/(rescale^(1/beta_v));
    E_zm0 = E_zm0/(rescale^(1/beta_m));        
    E_zm1 = E_zm1/(rescale^(1/beta_m));      

    
%% Solve the equilibrium    

    % Setting
    tol = 1e-12;
    maxIter = 1000;

    % Initial guess
    D_Hnew = ones(1,Q); 
    c_new = ones(1,Q);

    % Iteration
    err = 1;
    iter = 0;
    while err >= tol && iter < maxIter

        % new guess
        D_H = D_Hnew;
        c = c_new;

        % D(q,E)
        D0 = D_H;
        rv_ads = (D_H/f_vH).^(1/(beta_v-1))./((D_H/f_vH).^(1/(beta_v-1))+(e^sigma*D_F/f_vF).^(1/(beta_v-1)));
        D1 = rv_ads.*D_H + (1-rv_ads)*e^sigma.*D_F;        
        
        % Pi(q)
        Pi0 = D0.^gamma.*(c.^alpha_m*P_s^alpha_s).^((1-sigma)*gamma).*C_x0;
        power = beta_v/(beta_v-1);
        DE = (D_H.^power+(e^sigma*D_F).^power).^(1/power);
        Pi1 = Pi0.*(DE./D_H).^gamma;      
        
        % V(q) and M(q)
        f_v1 = rv_ads.^beta_v*f_vH + (1-rv_ads).^(beta_v)*f_vF;
        C_v1 = (1./(sigma*f_v1*w)).^(1/beta_v);       
        Vbar = C_v0*Pi0.^(1/beta_v).*E_zv0+rv_ads.*C_v1.*Pi1.^(1/beta_v).*E_zv1;      
        V = sum(phi_v.*repmat(Vbar,Q,1),2)';
        M = C_m*(Pi0.^(1/beta_m).*E_zm0+Pi1.^(1/beta_m).*E_zm1);

        % theta_v(q) and theta_m(q)
        xi = M./V;
        theta_v = 1-exp(-kappa*xi);
        theta_m = (1-exp(-kappa*xi))./xi;
        theta_m(isnan(theta_m))=kappa;       
        theta_m(theta_m==0)=kappa;
        
        % X(q,E)
        X0 = Pi0.*E_zx0;
        X1 = Pi1.*E_zx1;
        X_m = alpha_m*(sigma-1)/sigma*(X0+X1);

        % P(q)^(1-sigma)
        Psig = (X0./D0+rv_ads.*X1./D1);
        %Psig2 = X0./D0 + D_H.^(1/(beta_v-1))./(D_H.^(beta_v/(beta_v-1))+(e^sigma*D_F).^(beta_v/(beta_v-1))).*X1;
        %sum(abs(Psig - Psig2))
        
        % D_m(q)
        TempD = theta_v./M.*c.^(sigma-1).*X_m;
        TempD(isnan(TempD))=0;
        D_m = sum(phi_y.*phi_v.*repmat((TempD)',1,Q),1);

        % X_s
        rv_rev = rv_ads.*D_H./D1;
        %rv_rev2 = D_H.*D_H.^(1/(beta_v-1))./(D_H.^(beta_v/(beta_v-1))+(e^sigma*D_F).^(beta_v/(beta_v-1)));
        %sum(abs(rv_rev - rv_rev2))
        B1 = min(sum((1-rv_rev).*X1), 0.7); %(1-(sigma-1)*alpha_m/sigma) = 0.7360
        X_s = (1-(sigma-1)*alpha_m/sigma)*X-B1; 
               
        % D_s(q)
        D_s = phi_s.*X_s/sum(phi_s.*Psig);

        % D(q) (new guess)
        D_Hnew = D_m+D_s;

        % c(q) (new guess)
        c_new = (theta_m./V.*sum(phi_y.*phi_v.*repmat(Psig,Q,1),2)').^(1/(1-sigma));   
        
        % update
        D_Hnew = D_H + 0.2*(D_Hnew - D_H);
        c_new = c + 0.2*(c_new - c);
        iter = iter+1;      

        % check convergence
        err1 = max(abs(D_Hnew-D_H)./D_H);
        err2 = max(abs(c_new-c)./c);
        err = err1 + err2;

    end
    
    % Warning
    if iter == maxIter
        msg = 'User Warning: Inner Loop Maximum Iterations Reached.';
        warning(msg)
    end 

    % Adjust back
    Pi0 = Pi0/rescale;
    Pi1 = Pi1/rescale; 
    
    % Save the Export Decision
    ExportCutoffQ = (k0*lnzQ-log(sigma*gamma)+log(kron(ones(Nf,1),Pi1)-kron(ones(Nf,1),Pi0))-mu_E)/sigma_E;
    ExportProbQ = normcdf(ExportCutoffQ);        

    % Save mbar (in service price)
    P_sF = P_s/e*(B1/(X_s+B1))^(1/(1-sigma_s));    
    mbar = sum(V)/sum(phi_s.*Psig)*(P_s^(1-sigma_s)-(e*P_sF)^(1-sigma_s))^((1-sigma)/(1-sigma_s));
    f = X/(sigma*gamma*P_s);
    
%% Return

    eqm.Pi0 = Pi0;
    eqm.Pi1 = Pi1;
    eqm.D_H = D_H;
    eqm.D_F = D_F;
    eqm.D0 = D0;
    eqm.D1 = D1;
    eqm.D_s = D_s;
    eqm.D_m = D_m;
    eqm.c = c;
    eqm.Psig = Psig;
    eqm.P = Psig.^(1/(1-sigma));
    eqm.X_m = X_m;
    eqm.X0 = X0;
    eqm.X1 = X1;
    eqm.V = V;
    eqm.Vbar = Vbar;
    eqm.M = M;
    eqm.xi = xi;
    eqm.theta_v = theta_v;
    eqm.theta_m = theta_m;  
    eqm.rv_ads = rv_ads;
    eqm.rv_rev = rv_rev;
    eqm.C_v1 = C_v1;    
    eqm.B1 = B1;
    eqm.P_s = P_s;
    eqm.e = e;
    eqm.X_s = X_s;
    eqm.X = sum(X0+X1);
    eqm.P_sF = P_sF;
    eqm.mbar = mbar;  
    eqm.f = f;
    eqm.ExportCutoffQ = ExportCutoffQ;
    eqm.ExportProbQ = ExportProbQ;
    eqm.iter_inner = iter;
    eqm.err_inner = err;
    
    F = eqm;

    
end